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We numerically study the effects of two forms of quenched disorder on the anyons of the toric 
code. Firstly, a new class of codes based on random lattices of stabilizer operators is presented, and 
shown to be superior to the standard square lattice toric code for certain forms of biased noise. It 
is further argued that these codes are close to optimal, in that they tightly reach the upper bound 
of error thresholds beyond which no correctable CSS codes can exist. Additionally, we study the 
classical motion of anyons in toric codes with randomly distributed onsite potentials. In the presence 
of repulsive long-range interaction between the anyons, a surprising increase with disorder strength 
of the lifetime of encoded states is reported and explained by an entirely incoherent mechanism. 
Finally, the coherent transport of the anyons in the presence of both forms of disorder is investigated, 
and a significant suppression of the anyon motion is found. 
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I. INTRODUCTION 

A working quantum computer performing meaningful 
calculations unarguably requires information processing 
to be carried out in a fault-tolerant manner [TJ [2] . This 
not only means protecting the information from the ac- 
tion of imperfect gates, but also storing it in a reliable 
way during the course of computation. In the theory 
of quantum error correction, the state of a logical qubit 
can be encoded in the code space of a number of phys- 
ical qubits [3]. The resulting redundancy allows one to 
implement fault-tolerant quantum gates and to periodi- 
cally check for the occurrence of single-qubit errors using 
syndrome measurements. However, this kind of active 
error monitoring imposes an additional overhead on an 
already deeply involved vision. Therefore, the idea to 
manipulate and store quantum states in systems that al- 
ready provide 'built-in' protection from errors has gained 
a lot of attention recently [3HZ] ■ A promising approach in 
this direction is to encode information in the degenerate 
topologically ordered ground state of a suitable many- 
body Hamiltonian. Information is encoded in an entan- 
gled state distributed across a large number of qubits and 
can only be distinguished and modified non-locally. In 
this context, Kitaev's toric code [4j is arguably the best 
studied model to date. It is robust against local pertur- 
bation at zero temperature, as well as against thermal 
errors if long-range interaction between its fundamental 
excitations is present [8l |9] . 

Recent studies have focussed on coherent phenomena 
in the toric code that arise due to the additional pres- 
ence of various forms of quenched disorder |10H13j . Con- 
versely, this work is primarily concerned with a numerical 
study of incoherent (classical) effects caused by two par- 
ticular forms of randomness. First, we consider a class 
of models similar to the toric code, but differing from 
the latter in that the syndrome check operators are cho- 
sen randomly from a set of 3-body and 6-body opera- 
tors. These correspond to a generalization of the (square 



lattice) toric code to randomized lattices. We find that 
these models have an advantage over the toric code for 
biased noise, where bit-flip and phase-flip errors occur 
with different probabilities. We also present strong evi- 
dence that these codes are almost optimal, in the sense 
that they reach error thresholds close to the overall up- 
per bound valid for any Calderbank-Shore-Steane (CSS) 
code [Tallin]. Second, we investigate the effect of ran- 
dom onsite potentials on the lifetime of states encoded 
in the regular toric code coupled to a thermal bath. We 
identify and describe an interesting regime, where, in 
the presence of long-range interactions, the lifetime of 
this quantum memory is enhanced for increasing disor- 
der strength. Finally the effects of the random lattices 
on coherent anyon transport is investigated, both with 
and without additional randomness in the onsite poten- 
tials. The resultant slowdown of the anyonic motion is 
determined and its effect on the stability of the quantum 
memory is discussed. 

The paper is organized as follows: Section |ll] reviews 
the toric code, which is the basis of all further studies in 
this work. We then show in Sec. Illll how to simulate the 
classical dynamics of excitations in the systems consid- 
ered subsequently and also give some details on the nu- 
merics. Our main results for the random lattice models 
and random onsite potentials are presented in Sections 
IV and [V| respectively. Then, Sec. VI studies coherent 
anyon transport on the random lattices. A conclusion of 
our work is given in Sec. |VII| 



II. REVIEW OF THE 2D TORIC CODE 

The starting point of our investigation is Kitaev's 2D 
toric code [4 which will be modified in the following sec- 
tions to incorporate randomness. We provide here a brief 
outline of the original model for the sake of completeness. 
The 2D toric code consists of 2L^ spins- ^ with each spin 
placed on an edge of an underlying L x L square lattice 
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with periodic boundary conditions. One then defines two 
sets of mutuaUy commuting four-body operators, caUed 
plaquettes and stars, respectively, in the following way 
(cf. Fig. [IJi). A plaquette is the product of the Pauli 
a 2. operators associated with the four spins belonging to 
a single face of the square lattice, whereas a star is the 
product of the four cr^; operators of the spins on edges 
adjacent to a single vertex of the lattice. In this way, one 
obtains two sets of I? plaquette and star operators, out 
of which — 1 in each set are independent. Note that 
these operators can only have eigenvalues +1 and —1. 

One can then define a subspace C of the total Hilbert 
space given by the "^■'1 = 4 states which are si- 

multaneous eigenstates of all independent plaquettes and 
stars with eigenvalue +1. This space can thus accommo- 
date two logical qubits, and measuring the plaquette and 
star operators allows one to gather information about 
possible spin- and phase-flip errors without disturbing 
the encoded state. A negative plaquette (star) indicates 
the presence of one or three Ox is^z) spin errors. The toric 
code belongs to the class of stabilizer codes, and the pla- 
quettes and stars are in that context often referred to 
simply as stabilizer operators. 

Notably, the code space C is the degenerate ground 
space of the Hamiltonian 
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Here, J > is the energy gap, and A 
stars and plaquettes, respectively, explicitly given by 
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where adj(s) [adj(p)] denotes the set of spins on edges 
adjacent to the star (plaquette) s (p). The elementary 
excitations of the Hamiltonian Eq. ([T]) are stabilizer op- 
erators with negative eigenvalue and are referred to as 
'anyons'. 

Associated with the two logical qubits encoded in C are 
four string- like operators (products of single-spin Ox 's or 
cr^'s) which wrap around the torus and commute with all 
plaquettes and stars, but act non-trivially in the form of 
logical Pauli X and Z operators on the two qubits en- 
coded in C. We choose to label the operators such that 
Xx is a vertical string on horizontal edges and is a hor- 
izontal string on horizontal edges. Correspondingly, X2 
and Zi are horizontal and vertical strings, respectively, 
on vertical edges (see Fig. [l^). 

When the system described by the Hamiltonian Eq. ([T]) 
is coupled to a noisy environment causing single-spin er- 
rors, pairs of anyons are created and can subsequently 
move diffusively on the toric surface. Eventually, the cre- 
ation and diffusion of anyons leads to a pattern of errors 
containing undetectable loops around the torus, acting as 
unnoticed logical Pauli operators on the code space C and 
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FIG. 1. (color online) Toric codes, a) Kitaev's original 2D 
toric code. Shown is a 4 x 4 subregion of the L x L lattice. 
The blue solid dots on the edges of the lattice represent spins, 
the 4-body plaquette and star operators are shown in light 
green and orange, respectively (note that stabilizers contain- 
ing spins outside the figure are not shown). The four thick 
horizontal and vertical lines represent the logical Pauli oper- 
ator strings Xi,2 (red) and Zi,2 (green), b) The same region 
after introducing a regular pattern of spin defects. The modi- 
fied plaquette and star operators are not shown yet, see Fig. [2] 
The empty circles indicate the defect positions, i.e., the edges 
of the lattice where spins are removed. This requires altering 
the logical operators Z2 and X2 in the way shown. Note that 
all commutation relations between the logical operators are 
preserved. 



therefore corrupting the state contained therein. Mea- 
suring the plaquette and star operators to locate anyons 
reveals some, but not all information about the under- 
lying error pattern and is generally ambiguous. It is up 
to an error correction procedure (see Sec. 
with this problem in a satisfactory way. 
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The toric code has gained attention due to a series 
of interesting and advantageous properties. Namely, the 
stabiUzer operators are local and independent of the sys- 
tem size L, while the code distance grows linearly with L. 
Closely related is the fact that the ground-state degen- 
eracy is exponentially protected (in L) against local per- 
turbations. Quite remarkably, the toric code is in a sense 
almost optimal within the class of all CSS codes even 
though the latter contains codes with arbitrarily large 
stabilizer operators. We will revisit this topic in greater 
detail in Sec. HVCl 



time scale, /3 — l/fcsT, with T being the temperature of 
the bath and ks denoting Boltzmann's constant, and uJc 
is the cutoff frequency of the bath. In the following, we 
set UJc ^ oo for simplicity. A bath with rt = 1 is called 
'Ohmic', whereas one with 71 > 2 is called 'super-Ohmic'. 
Only the former case is considered in this work. Unless 
otherwise stated, all energies will be expressed in units 
of ksT. Consequently, the unit of time is (ki/csT)"^. 



B. Simulations and error correction 



III. CLASSICAL DYNAMICS AND 
NUMERICAL SIMULATIONS 

A. Classical dynamics from single-spin errors 

Since the Hamiltonian Eq. ([T]) does not couple the 
star and plaquette operators, we can treat the two cor- 
responding types of anyonic excitations independently. 
Furthermore, because the stars are simply plaquettes on 
the dual lattice, it is sufficient to study the dynamics of 
only one type, e.g., the plaquette anyons. We assume 
that each spin is coupled to an auxiliary system that can 
cause the spin to flip via ax errors. In the limit of weak 
coupling [THl ET] , one can derive the following system of 
coupled rate equations describing the classical dynamics 
of the system |^ : 
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Here, pe{t) is the time-dependent probability to find the 
system in the state \£) obtained by applying errors to 
all spins with indices in £, i.e., \£) = Yikee '^x\'^o)i where 
\^o) is the initial state of the system. Similarly, Pxi(e)i^) 
describes the probability to be in the state cr^|f ). Finally, 
7* £ and 7°£* are the transition rates to arrive at or leave 
the state \£), respectively, via a CT2;-error at the spin with 
index i. 

In this work, we will consider two types of error envi- 
ronments. The first one is a constant error rate model, 
i.e., we set 7*'^ = 7™* = const. In this case, spin- 
flips are caused independently of any previously exist- 
ing anyons and Gx errors. The second model mimics 
the coupling to a thermal environment, where the tran- 
sition rates are in general energy dependent. Conse- 
quently, we set 7^'"^ = -f{-uji^s) and 7°J* = ^{uJi^g), where 
^i,S = ^£ ~ ^x^(£) is the energy difference between the 
states \£) and cr^lf). An expression for 7(aj) often found 
in the literature is given by 



7(w) = 2k„ 
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and can be derived from a spin-boson model [181 I19j . 
Here, constant with units 1/energy" setting the 



Clearly, it is impossible to solve Eq. Q analytically for 
meaningful system sizes, because the number of states pg 
grows exponentially with L^. We thus have to stochas- 
tically simulate the system and obtain the quantities of 
interest, such as the number of anyons or the expecta- 
tion values of the logical operators, by averaging over 
many (typically several thousand) instances. In greater 
detail, the iteration of the simulation at time t consists of 
these steps: (i) Calculate all unnormalized single spin-flip 
probabilities pi — j{e£ — £xi{£))i then obtain from them 
the total spin flip rate R = 'Y^^Pi- (ii) Draw the time 
At until the next spin flip from an exponential distribu- 
tion with rate R. (iii) Calculate and record all quantities 
of interest for time sampling points lying in the interval 
[t,t + At]. Namely, these quantities are the number of 
anyons, the number of Ux errors, and the uncorrected and 
error-corrected (see below) logical operators Zi and Z2. 
(iv) Determine a random spin according to the probabil- 
ities Pi/R, flip it, and set t to t + At. 

The error correction step in the toric code consists 
of pairing up all detected anyons and then annihilating 
them by connecting each pair with a string of errors from 
one partner to the other. The pairing is usually chosen 
such that all anyons are annihilated with the smallest 
total number of single-spin operators. This is known as 
the minimal-weight perfect matching and can be found 
in polynomial time with the help of the 'blossom' algo- 
rithm due to Edmonds |20j . The runtime complexity 
of this algorithm has been improved several times since 
its discovery. We are employing the library Blossom 
V 21J which implements the latest version running in 
0{mnlogn) time, where n is the number of anyons (ver- 
tices) and m the number of connections (edges) between 
them. 

In order to find the true matching with minimal weight, 
one in principal would need to choose the set of edges to 
include all connections from every anyon to every other. 
However, since the size of this set grows quadratically 
with the number of anyons n, the overall scaling of the 
matching algorithm becomes 0{n^ \ogn) which is infea- 
sible for large n. We therefore first perform a Delaunay 
triangulation in negligible O(nlogn) time using the li- 
brary Triangle [22] . The result is that only anyons close 
to each other are connected using a number of edges lin- 
ear in the number of anyons. It turns out that this is an 
excellent approximation yielding results that are nearly 



indistinguishable from those obtained from a matching 
over the complete graph. 

Within the paradigm of active error correction, where 
the anyons are detected and corrected periodically on 
sufhciently small time intervals, the encoded state can 
be kept free of logical errors almost indefinitely. How- 
ever, since we are interested in the use of the toric code 
as a passive quantum memory, we are mostly concerned 
with the lifetime r of the encoded information in a sce- 
nario where error correction is only performed once at 
readout. Hence, whenever we show plots of the 'error- 
corrected' logical operators decaying as a function of 
time, we thereby refer to their value if error correction 
had been performed at that time, without actually per- 
forming it. We then define the lifetime of the system as 
the time it takes for the expectation values of the error- 
corrected logical operators to decay to 90% of their initial 
value. 



IV. RANDOM LATTICES 

In this section, we study the error thresholds of a fam- 
ily of models obtained by randomly modifying the toric 
code in a way that preserves its basic features. We first 
describe how we create our random lattices and then 
present and discuss the results of the simulations within 
the context of optimal quantum codes. 




A. Generating the lattices 

Starting from the toric code on an L x L lattice, we 
remove ^L^ spins at specific and regularly distributed 
'defect' locations. The structure of the defect pattern 
can be easily understood from Fig [TJd. Basically, every 
second vertical edge is labelled as a defect, with the first 
defect of each row alternatingly being created on the first 
or second vertical edge of that row. Note that the height 
and width of the grid both must be even in oder for this 
procedure to be consistent with the periodic boundary 
conditions. We now have to modify all plaquettes and 
stars as well as the logical Pauli X and Z operators such 
that all original commutation relations remain unaltered. 

Let us start with the logical operators. Clearly, both 
Xi and Zi (with single-spin operators exclusively on hor- 
izontal edges) are unaffected by the introduction of de- 
fects on vertical edges of the lattice. However, the pair 
of operators in the original code acting on the second en- 
coded qubit is defined on vertical edges and thus needs 
to be adapted. Clearly, the operators must remain con- 
nected strings wrapping around both dimensions of the 
torus. The 'zig-zag' pattern shown in Fig. [T|d achieves 
this with the smallest increase in the number of single- 
spin operators. It is straightforward to verify that these 
new X2 and Z2 operators, together with the unaltered 
pair acting on the first encoded qubit, indeed fulfill all 
original commutation relations between each other. 



FIG. 2. (color online) Modifying the stabilizer operators, a) 
When removing a spin (empty circle), we choose between two 
ways of adapting the affected stabilizer operators. With prob- 
ability Pmix we join the two plaquette operators to one large 
6-body operator and reduce the two stars from 4- to 3-body 
operators. Alternatively, with probability 1 — Pmix, we per- 
form the dual operation, namely we define two 3-body pla- 
quettes and one large 6-body star. Spins and operators not 
affected by removing the central spin are not shown in this 
example, b) Typical 8x8 subregion of a (larger) random lat- 
tice with Pmix = 0.5. Logical operators as well as some spins 
and operators at the edges with the rest of the lattice are not 
shown. 



We now discuss the modification of the plaquette and 
star operators. Removing one spin, i.e., creating one de- 
fect, affects exactly two adjacent plaquettes and two ad- 
jacent stars. We will consider two possible ways of deal- 
ing with this situation (cf. Fig. [2]). We can either (i) 
define two restricted 3-body plaquette operators and one 
'vertical star' operator consisting of the product of the 
remaining 6 single-spin ax operators, or (ii), perform the 
dual operation, namely defining one large 6-body 'hor- 
izontal plaquette' and two restricted 3-body stars. The 
two ways of modifying the original operators are depicted 
in Fig [2^. It is relatively easy to see that these new 3- 
body and 6-body operators remain mutually commuting, 
and furthermore commute with all modified logical Pauli 
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operators just as in the original code. Note that also the 
dimension of the code space is left unchanged since it 
generally only depends on the genus of the surface cov- 
ered by the anyon operators We can now create a 
random lattice by choosing at each defect site to cre- 
ate a 6-body plaquette with probability Pmix and two 
3-body plaquettes with probability 1 — Pmix ■ See Fig. [SJd 
for a typical example. The special case Pmix = cor- 
responds to a regular lattice of 3-body plaquettes and 
6-body stars, whereas Pmix = 1 conversely yields a reg- 
ular lattice of 6-body plaquettes and 3-body stars. In 
both cases, the 3-body operators are the vertices of an 
underlying hexagonal lattice, whereas the 6-body opera- 
tors form the vertices of its dual, the triangular lattice. 



B. Results 

The randomization of the lattice geometry described 
above creates a difficulty in employing the Delauny tri- 
angulation, because the graph is now irregular. To 
overcome this problem we have replaced this step by a 
breadth-first search performed on each anyon. This pro- 
cedure connects every anyon to at most k of its nearest 
neighbors, where distance is measured not in a Euclidean 
sense, but as the number of errors in a connecting string. 
For constant k, this requires a runtime of 0(n), where n 
is the number of anyons. We have found that fc = 10 is 
an excellent approximation to fc = cx) and have used this 
value in all calculations. 

We have performed a series of Monte Carlo simulations 
to determine the critical fraction of errors /J! indepen- 
dent of any form of anyon dynamics (see Appendix [A| for 
additional results in the case of thermal errors). If the 
probability for each spin to independently be affected by 
a ax error becomes larger than this critical value in the 
limit L — >■ oo, the error correction scheme undergoes a 
transition from performing fully accurate error recovery 
to randomly guessing the error-corrected state with the 
lowest possible success rate of 50%. We can determine 
by plotting the expectation values of the error-corrected 
logical Z operators as a function of the error probability 
/ for different lattice sizes and observe at which value of 
/ the curves intersect [25] . 

Fig. [3^ displays typical results for a few different val- 
ues of Pmix- We observe that the expectation values of 
Zi and Z2 in general have a different dependence on the 
error probability /. This is due to the fact that the frac- 
tion of 6-body plaquettes determines the average number 
of spins required to form a loop around the horizontal 
direction of the torus, whereas this number is constant 
for the vertical direction. Not surprisingly, Zi and Z2 
decay identically for a 50 per cent mixing of 3- and 6- 
body plaquettes. Note that, despite the typically un- 
equal decay of Zi and Z2 as a function of /, the error 
thresholds for the two operators are always identical for 
a given value of Pmix- This is consistent with the gen- 
eral understanding that the correctability of the memory 
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FIG. 3. (color online) Critical error thresholds of random 
models, a) Three example plots of data used to determine 
the critical error thresholds. Each plot shows, for a spe- 
cific value of Pmix, the expectation values of the logical Zi 
(dotted lines) and Z2 (dashed lines) operators for grid sizes 
L — 32 (circles), 64 (triangles), and 128 (squares). The ver- 
tical dotted lines indicate the position of the error threshold. 
Data points are obtained by bootstrapping 1000 sample val- 
ues, each of which is obtained by averaging over 200 random 
error distributions on a single instance of a random lattice, 
b) Error thresholds [as determined in a)] of Z (circles) and X 
(triangles) operators as a function of pmix ■ The dotted lines 
are guides to the eye. 
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as a whole is related to the phase of a corresponding 
random-bond Ising model [5 . Indeed, our numerically 
determined thresholds for the regular 3-body and 6-body 
plaquette lattices agree well with the recently calculated 
multicritical points in spin glass models on hexagonal 
and triangular lattices, respectively [H]. Specifically, 
we find fc ~ 0.1585 for Pmix = (theoretical value: 
fc = 0.1640) and fc ~ 0.0645 for Pmix = 1 (theoretical 
value: fc = 0.0674). The discrepancy of about 4 — 5% 
between the numerical and theoretical thresholds is of 
the same size as in the case of the toric code on a square 
lattice (where we had found fc « 0.1055 as compared to 
the theoretical value fc — 0.1092) [9^ and can generically 
be attributed to the failure of the minimal-weight per- 
fect matching close to the threshold. Interestingly, our 
numerical results suggest that the thresholds of the toric 
code and our random lattice models with Pmix — 0-5 are 
identical. 

We show in Fig.jsjD the critical fraction of errors /J! for 
the logical Z operators determined in the way described 
above as a function of the lattice mixing probability Pmix ■ 
Since the plaquette and star lattices are dual to each 
other (to every 6-body plaquette correspond two 3-body 
stars and vice versa) , the critical fraction f^ of Cz errors 
for which error correction of the logical X operators fails 
(also plotted in Fig. ^) is simply given by 



(6) 



fcriPniix) — fcri^ Pmix)- 



At equal mixing, i.e., Pmix — 0-5, the threshold values 
are given by /^f (0.5) = ^^(0.5) « 0.1055. 

Consequently, one of the thresholds for the two differ- 
ent types of Pauli operators, either or /J! is always 
smaller than or equal to the threshold of the toric code. 
Our random lattices thus bear no advantage over the lat- 
ter in the case of a uniform error model, where cr^ and 
Gz errors occur with the same probability. The situa- 
tion is different, however, for biased noise. If bit-flips 
and phase-flips are created with different probabilities, 
we can make use of the asymmetry in the error thresholds 
ioT Pmix 7^ 0-5- Assuming, for instance, that errors are 
more frequent than Uz errors would lead to an overall life- 
time decrease of encoded states in the toric code due to 
the shorter lifetimes of the logical Z operators. However, 
starting from a random model at Pmix = 0.5, a decrease 
in Pmix will lead to an increase of the Z lifetimes and a 
decrease of the X lifetimes. If the error frequencies are 
not too different, the lifetimes will become identical at 
some value < Pmix < 0.5 and will be larger than the 
overall lifetime of the toric code. We thus conclude that 
the random lattices can be employed to increase the life- 
time of encoded states compared with the toric code on 
a square lattice in the presence of biased noise. While 
these lattices require both error probabilities to fall in 
the range 0.0674 < p < 0.1640 (and below the bound- 
ary in Fig. |4j see next section) , it should in principle be 
possible to extend this range to < p < 0.5 by defining 
stabilizers with more than 6 single-spin operators in a 
similar fashion. 




FIG. 4. (color online) Theoretical upper bound on biased 
noise correctable by CSS codes. Given an error model of 
independent and az errors occurring with constant proba- 
bilities Px and Pz, respectively, there exists no CSS code able 
to cope with pairs of error probabilities lying above the solid 
line given by the zero-contour of Eq. ([7|. The crosses are 
the numerically determined pairs of thresholds of the random 
models for Pmix ~ 1 down to Pmix = from left to right. The 
dotted line is a guide to the eye. 



C. Relation to optimal quantum codes 

We now discuss our random lattices within the con- 
text of optimal quantum coding. It is well known that, 
assuming a biased constant error model, there is an up- 
per bound on the fraction of logical qubits k and physical 
qubits n that encode them, valid for all CSS codes. This 
bound is given by [51 [151 Hi] 



fc/n < 1 - F(p,) - Hipz), 



(7) 



where H{x) — —x log2 x~(\ — x) log2(l — x) is the Shan- 
non entropy, and px and Pz are the probabilities for a 
single spin to be affected by a and Oz error, respec- 
tively. The bound Eq. ^ can be motivated with the 
following intuitive argument. 

An ideal CSS code would be able to detect for each 
physical qubit if it was suffering from a Gx or a Gz error. 
Assuming that these errors are uncorrelated, the num- 
ber of classical bits required to store this information is 
asymptotically given by nH(px) + nH{pz). If we are to 
store the same information in qubits instead of bits, the 
Holevo bound [1 requires the usage of at least as many 
qubits to do so. Since our optimal CSS code needs to 
store the information of k encoded qubits, as well as all 
possible occurrences of errors, we have 



> k + nH{px)+nH(pz). 



(8) 



Dividing by n and rearranging the terms yields the de- 
sired bound Eq. ([7]). 

For unbiased errors with p = Px ~ Pz, the right-hand 
side of Eq. ([t]) becomes zero at p w 0.110028, implying 
that there cannot be any CSS code coping with an er- 
ror rate larger than this value. Quite remarkably, the 
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critical error probability for the toric code has been de- 
terniined to be fcr — 0.109187 [23]. This is astonishingly 
close to the upper bound, especially when taking into 
account that all stabilizers are local 4-body operators. 
Moreover, inserting the error thresholds for the regular 
6-body plaquette lattice and its dual lattice of 3-body 
stars {pmix = 1) into the right-hand side of Eq. ([t]) eval- 
uates to 

1 - iy(0.0674) - iJ(0.1640) w 6 x 10-^ (9) 

which is virtually zero, indicating that the code is close to 
optimal for this particular biased error model. Due to the 
symmetry of Eq. Q with respect to the error probabili- 
ties and the duality of the triangular and hexagonal lat- 
tices, the same argument holds for a lattice with 3-body 
plaquettes and 6-body stars {pmix = 0) with the values 
of Px and Pz exchanged. With our random lattices, we 
can thus continuously interpolate between two optimal 
models by changing Pmix ■ This suggests that the random 
models are optimal for all values of Pmix , in the sense that 
for every 0.0674 < Px < 0.1640 there is a random model 
with a theoretically (close to) maximal possible thresh- 
old for Pz ■ The results plotted in Fig. [4] strongly support 
this claim. The solid line is the zero-contour of the upper 
bound Eq. ([T]) and the crosses are the threshold pairs of 
the random lattices determined numerically. Note that 
the numerical data is within the typical 5% distance of 
theoretical bound. This can once again be explained by 
the failure of the minimal-weight error correction algo- 
rithm close to the thresholds. This observation, together 
with the knowledge from theory that the models are vir- 
tually optimal for Pmix = and Pmix = 1 leads us to 
conjecture that the random models are virtually optimal 
for all values oipmix- However, carrying out a theoretical 
study in the fashion of Ref. ^24] is outside of the scope of 
the present work and is deferred for future research. 



V. RANDOM ONSITE POTENTIALS 

This section is devoted to the study of the classical dy- 
namics of anyons in the regular toric code on a square lat- 
tice, but with randomly modified anyon onsite energies. 
We are also particularly interested in the case where long- 
range anyon-anyon interaction is present, as this has been 
shown to generally enhance the lifetime of the memory 
due to the suppression of the anyon density with increas- 
ing system size [9j. For this, it is convenient to intro- 
duce the new stabilizer operators rig = (1 — Ag)/2 and 
Hp = (1 — Bp)/2, where and Bp are the usual star and 
plaquette operators, respectively. These operators are 
zero in the absence of an anyon on the respective site, 
and equal to 1 otherwise. The more general Hamiltonian 
can then be written as 

H ^^^^^Upp^UpUp^ + ]^^^Vss'nsns', (10) 

pp' ss' 



where Uppi and Vss' contain the onsite energy and repul- 
sive anyon interaction terms. Since in this model pla- 
quette and star anyons are still independent, we can set 
Vss' = and note again that all results for the plaquette 
anyons hold equally for the stars. We then set [9] 

Upp' = 2Jp6pp/ -f - ~ Spp')i (11) 

Vpp' ) 

where Jp is the onsite energy of an anyon on the plaque- 
tte with index p, A is the interaction strength, rppi is the 
shortest distance on the torus between plaquettes p and 
p' , and < a < 2 is the (long-range) interaction expo- 
nent. The onsite energies Jp are chosen randomly from 
a distribution with mean zero in order to discriminate 
effects caused by the randomness from effects potentially 
caused by the system having a non-zero mean gap. 

We focus on the case of constant interaction, i.e., 
a = 0, and Ising-like randomness, meaning that the Jp 
are chosen from {— <t, +(t} with equal probabilities. We 
refer to tr > as the disorder strength. This model is in- 
teresting mostly for two reasons. First, it is the most con- 
venient system incorporating randomness with respect to 
numerical simulation. Since the interaction is constant 
and thus simply depends on the total number of anyons, 
only six different single-spin flip rates need to be updated 
at each iteration step, depending on the number and con- 
flguration of adjacent anyons and onsite energies, respec- 
tively. Second, this simple model already displays all 
dynamical effects also present in more complicated sys- 
tems (e.g., a ^ and Gaussian distribution of Jp's, see 
Appendix |B]) and thus serves as an ideal playground for 
studying these effects. Naively, one would expect that 
the presence of negative onsite energies in the system 
simply favors the creation of anyons and is thus always 
disadvantageous for the lifetime of the memory. While 
this is indeed true for a non-interacting system, we find a 
regime in the interacting case where, quite surprisingly, 
the lifetime is enhanced for increasing disorder strength. 

Fig. |5] presents the results for the two cases. The non- 
interacting system is stable against disorder strengths 
that are roughly equal to the temperature but then de- 
cays for larger a. This can be understood easily from the 
detailed balance condition satisfied by the rates Eq. ([5|: 
The ratio of the creation and annihilation rates of a pair 
of anyons on two sites with negative onsite energy is given 
by 7(2<t)/7(— 2(7) = cxp(2/?cr), which becomes large for 
ct's exceeding the temperature. It is thus exponentially 
more likely for a pair of anyons to be created than annihi- 
lated for (T > fc^T, thereby quickly cluttering the system 
with anyons and crossing the critical fraction of errors. 
This situation changes completely in the presence of in- 
teractions between the anyons. The data shown in Fig. [5] 
displays a steep increase in the lifetime as a function of 
the disorder strength, peaking at around a « Z.bksT for 
that particular system, followed by a slower decay. 

We can shed some light on this effect by additionally 
looking at the number of anyons and the number of er- 
rors, as shown in Fig. |6] For any fixed value of ti, the 
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FIG. 5. (color online) Influence of disorder on the memory 
lifetime, a) Time evolution of error corrected Z operators for 
a non-interacting (top) and interacting (bottom, a = 0, A = 

0. 5) model. Both systems are of size 32 x 32 unit cells (2 x 32^ 
spins) and are coupled to an Ohmic bath at temperature T = 

1. The different curves display (Z^dt)) for different disorder 
strengths a of Ising-like randomness with onsite energies Jp = 
±a. In the non-interacting system, a is increased from 
to 10 ksT as indicated in the panel. The inset displays the 
lifetime of the memory, i.e., the time at which (Zee) hits 0.9, 
as a function of a. The disorder strengths examined in the 
interacting case have been chosen as < o" < 15fcsT (see 
main text and labels in the panel), b) The lifetimes of the 
interacting model extracted from the curves of the lower panel 
in a). The dotted line is a guide to the eye. 
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FIG. 6. (color online) Number of anyons (top) and number 
of single-spin errors (bottom) as a function of time in an in- 
teracting system {a — 0, A = 0.5 ksT) of size 32 x 32 coupled 
to an Ohmic bath. The strength of the Ising-like randomness 
is increased from cr = to cr = 15 fesT as indicated in the 
panels (see also main text). 



number of anyons again increases quickly but then sat- 
urates almost instantaneously at the equilibrium value. 
At this point, creating a new pair of anyons costs an en- 
ergy penalty due to the repulsive interaction that can 
no longer be compensated by the negative onsite poten- 
tials. One can clearly see that the equilibrium number 
of anyons increases linearly with a, which implies that 
the corresponding enhanced memory lifetimes cannot be 
explained by a suppression of the anyon density. 

However, the error creation rate (i.e., the slope of the 
curves in the lower panel of Fig. [6| exhibits a pronounced 
minimum at the same value a ~ i.bksT that also yields 
the maximal lifetime. Such an initial decrease in the er- 
ror rate despite an increasing number of anyons can only 
be consistently explained by a suppression of the anyon 
diffusion. For disorder strengths a 3> fcsT, processes 
that create anyons on two positive sites or that move 
an anyon from a negative to a positive site are exponen- 
tially suppressed. The positive sites thus effectively act 
as infinite barriers that greatly reduce the mobility of the 
anyons, and the encoded state is solely destroyed by dif- 
fusing anyons restricted to the negative sites [26] . As the 
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disorder strength is lowered, two competing effects come 
into play. On the one hand, the number of anyons is re- 
duced linearly. On its own, this would lead to a longer 
lifetime due to the presence of fewer diffusing anyons. 
On the other hand, the barriers separating the regions 
of negative sites are lowered, which facilitates the diffu- 
sion across longer distances and promotes a reduction in 
lifetime. The observed maximum in the memory lifetime 
can thus be understood as a tradeoff between having few 
but relatively freely moving anyons for a <^ kgT, and 
more but very restricted anyons for a 3> ksT. The in- 
teraction merely plays the role of restricting the anyons 
to a small enough (for a < ksT) and constant num- 
ber. Appendix [C] contains results that further support 
the picture described above. 

VI. QUANTUM DYNAMICS 

A. Error model 

Consider the toric code Hamiltonian, perturbed by a 
magnetic field of strength h. For concreteness, let us 
choose this to be of the form, 

H ^^JpUp + ^JsUs + h'^ai. (12) 

p s i 

The effects of such a perturbation have been studied us- 
ing the methods of Refs. [11] [T2], where it was noted 
that, since the al. do not commute with the rip, the per- 
turbation will have the effect of creating, annihilating 
and transporting plaquette anyons. For /i <C Jp all these 
effects apart from the transport are suppressed by the 
energy gap, allowing the system to be modelled as the 
following many-particle quantum walk Hamiltonian, 

Hp ^^Mpytpy + U ^UpiUp - 1), 

p,p' p 
Mp^p' = Si^ppi^h + Sp^p'Jp. (13) 

Here (5(p,p') = 1 only when the plaquettes p and p' share 
a spin. The operator tpy maps a state with an anyon 
on the plaquette p to that with the anyon moved to p', 
and annihilates any state without an anyon initially on p. 
Since the anyons are hardcore bosons, we are interested 
in the case of ?7 — >■ oo. 

This effective description in terms of quantum walks of 
anyons holds also for a more general magnetic field and 
other local perturbations. The effects of anyonic braiding 
occur at a higher order of perturbation theory than those 
of this effective description, and hence may be ignored. 
The dynamics of the plaquette and vertex anyons can 
therefore be considered separately. Since they are dual 
to each other, once again only the plaquette anyons are 
considered here without loss of generality. 

The Hamiltonian Hp is difficult to solve in general. 
However, note that the dynamics of Hp are driven by 
the matrix M, i.e., the Hamiltonian for a single particle 
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FIG. 7. Time evolution of the standard deviation A{t) of 
a single quantum walker on lattices with uniform couplings 
Jp = J for all p. The different curves correspond to a square 
lattice (solid), 3-body plaquette lattice (dashed), and 6-body 
plaquette lattice (dotted). Inset: A{t) on the random lattices 
from Sec. |IV| with Pmix = 0.5, where each point has been 
averaged over 1000 samples. 



walk. Hence, by considering the case of a single anyon, 
important aspects of the behavior for the many-particle 
walks can be determined. It is this approach that is taken 
here. The Hamiltonian M is applied to a single anyon, 
initially placed on an arbitrary plaquette of the code. 
The motion of the anyon can be characterized by the 
time evolution of its standard deviation, A. Since finite 
values of the system size L must be used in the numerics, 
the walks will, at some point, interact with the boundary. 
In order for this effect to be ignored, the run time of the 
walks is limited to ensure this interference always remains 
negligible. 

The behavior of A over time for walks on square, tri- 
angular and hexagonal lattices for which all Jp are uni- 
form can be found in Fig. [7j In each case the standard 
deviation of the distance increases linearly with time, 
demonstrating the ballistic motion expected from quan- 
tum walks when no disorder is present. 

The ballistic motion caused by the field is highly dam- 
aging to the quantum information stored within the code. 
Suppose that the toric code initially has some density p 
of anyon pairs, due perhaps to noisy preparation of the 
state or interaction with the environment. If p is suf- 
ficiently small then the pairs will be far apart, allowing 
error correction to be performed reliably. However, when 
the perturbation is present this will only remain true for a 
finite lifetime r, after which the motion of the anyons pre- 
vents them from being paired reliably. This occurs when 
they have moved a distance comparable to the average 
distance between pairs, and hence when A(t) ~ 
Since A(t) grows linearly with time for ordered quantum 
walks (cf. Fig. [7|, the quantum memory will fail within 
a time of order r = 0{l/y/p). Mechanisms which slow 
down the anyons are therefore favorable to the quantum 
memory, since they lead to longer lifetimes. It is this 
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effect that is expected to emerge from the disorder. 



B. Random lattices 

Let us now introduce disorder by using the random 
lattice of Sec. IV while still keeping the Jp (and Jg) uni- 
form, all taking the same value J. Specifically the case of 
Pmix — 0.5 is considered, to maintain the symmetry be- 
tween plaquette and star anyons. The behavior of the 
standard deviation of the distance moved by a single 
walker is shown for this lattice in the inset of Fig. [7] 
Rather than increasing linearly with time t, as in the or- 
dered case, it is found that A(t) grows with the square 
root of t. The motion of a quantum walker is therefore 
diffusive rather than ballistic in this case. As such, the 
random lattice leads to a significant slowing of the anyon 
motion, increasing the lifetime to r = 0[p~^) (note that 
we always have p < 1). It is possible that the random 
lattice also induces Anderson localization [57], in which 
case the lifetime will be increased further, but the system 
sizes which may be probed are too small for this to be 
evident. 
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FIG. 8. Time evolution of the standard deviation A{t) for a 
single anyonic walker with disorder in the Jp ot a/h = 250, on 
both a square lattice (solid line) and a random lattice (dashed 
line) with pmix ~ 0.5. Each point has been averaged over 1000 
samples. 



VII. SUMMARY 



C. Random lattices together with J disorder 

It is known that, when disorder in the Jp couplings 
is present in the toric code, Anderson localization is in- 
duced [TTJ [T^ . This effect exponentially suppresses the 
motion of the walkers, and causes the standard devia- 
tion of the distance to converge to a constant value. We 
therefore have t — )> cxd, i.e., the memory stays stable 
against the perturbation for an arbitrarily long time. It 
is now important to determine whether the combination 
of randomness in both the lattice and the Jp couplings 
enhances or diminishes this effect. 

To study this, disorder in the Jp couplings are con- 
sidered. Specifically, each Jp randomly takes either the 
value J — (7 or J + a with equal probabilities. The value 
J is unimportant, but the ratio of a/h characterizes the 
strength of the disorder in comparison to the magnetic 
field. Guided by the numerical results of [111, we con- 
sider here disorder of strength a/h = 250 to ensure that 
the localization effect is observed for moderately sized 
systems. 

In Fig. |8j the time evolution of the standard deviation 
is shown for the case of J disorder on a square and a ran- 
dom lattice. In both cases the localization effect is seen, 
with the walk unable to move far beyond a few times the 
length scale separating neighboring vertices. The walks 
with and without the lattice disorder give very similar 
results, especially at longer times. The effect of local- 
ization in the random lattice therefore seems the same 
as that of the square, without significantly enhancing or 
diminishing the effect. 



We have studied the infiuence of quenched disorder on 
the incoherent (classical) motion of anyons in modified 
forms of the toric code. We have first described a class 
of random models that can be obtained from the toric 
code by removing a regular sublattice of spins, and then 
for each defect site randomly choosing one of two ways 
to adapt the affected stabilizers with a probability Pmix ■ 
The critical fractions of independent errors at which these 
codes become uncorrectable have then been determined 
numerically as a function of Pmix- We have shown that 
in the presence of biased noise, where bit flips and phase 
flips occur at different probabilities, the models based on 
random lattices can tolerate higher thresholds than the 
toric code in one type of errors, given the other type 
is correspondingly lower. These thresholds have been 
demonstrated to be close to the upper bound correctable 
by any CSS code. Second, we have studied the toric code 
subject to randomness in the onsite potentials. Specifi- 
cally, we have demonstrated that in the presence of re- 
pulsive long-range interaction between anyons, there is a 
pronounced maximum in the lifetime of encoded states 
as a function of disorder strength. This effect has been 
attributed to a reduction of anyon diffusion due to the 
sites with positive onsite energy acting as barriers for the 
anyons. Finally, the effects of both forms of disorder are 
studied for coherent transport of the anyons. It is found 
that the random lattices cause the anyons to move diffu- 
sively rather than ballistically, increasing the lifetime of 
the memory. Adding randomness in the potentials then 
causes the anyons to localize, leading to further stability. 
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FIG. 9. (color online) Critical fraction of errors fcr ~ fit — r) 
as a function of lattice mixing Pmix at temperatures T — 0.2 J 
(circles), 0.45J (triangles), 1.0125J (squares), and 2J (dia- 
monds). The lifetime t is given by the intersection of error- 
corrected logical Z operators for lattice sizes L = 38, 56, 86. 
Error bars are due to the uncertainty in r. The solid line 
is determined by the thresholds from the Monte Carlo simu- 
lations of an independent error model (see main text). The 
inset shows an example of crossing logical Z operators for the 
particular values Pmix = 0.4 and T — 0.45J. 
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Appendix A: Critical fraction of random lattices in 
contact with an Ohmic bath 



We present here some additional results for the error 
thresholds of random lattices in the presence of thermal 
errors. Fig. |9] shows the fraction of errors / at the life- 
time of an infinitely large system as a function of the 
lattice mixing Pmix and for different temperatures T. In 
this section, the energy scale is set by the anyon gap 
J. Time is thus measured in units of (ki J)~^. For given 
Pmix and temperature T, we first simulate systems of sev- 
eral different sizes in contact with an Ohmic bath. We 
then determine the lifetime r as the intersection of the 
decay curves of the corresponding error-corrected logical 
Z operators (see inset of Fig. [9]). Since the anyon dy- 
namics are independent of the system size (note that the 
anyons are not interacting with each other), all curves 
f{t) for different system sizes collapse and the specific 
value f{t = r) = fcr can be read off easily. One can 
see nicely that these thresholds converge with increasing 
temperature to the ones given by the model of indepen- 
dent errors. This can be explained by the loss of cor- 
relations between errors due to an increasing amount of 
thermal noise in the form of fluctuating anyons. 



Appendix B: Gaussian noise and 1/r interaction 

We have also performed simulations with 1/r inter- 
action (a = 1) and plot the results in Fig. 10 Apart 
from Ising-like disorder (Jp = zttr) we have also looked 
at a Gaussian distribution of onsite potentials Jp with 
mean zero and standard deviation a. Generally, the life- 
times are shorter than for constant interaction because 
the weaker 1/r interaction allows for a larger density of 
anyons. Nevertheless, the results are qualitatively similar 
to the ones discussed in the main text, namely showing 
a pronounced maxima of the lifetime as a function of a. 
This supports the picture that the interaction is required 
to limit the number of diffusing anyons, while the ini- 
tial increase in lifetime with a is due to their obscured 
diffusion. In the case of Gaussian noise, the latter ef- 
fect is even stronger, because anyons are created or get 
trapped in a few sites with onsite potentials much lower 
than — (7, out of which it is difficult for them to escape 
again. This explains the increased lifetime from Ising to 
Gaussian randomness. 



0.75 



b- 0.5 



0.25 - 



1 1 1 1 1 1 1 1 


1 




'a 






A a 












' O . 


~~A. 




1 , 1 . 1 , 1 , 1 , 



2.5 5 7.5 10 12.5 15 

FIG. 10. (color online) Lifetime r as a function of disorder 
strength a in the presence of Gaussian (triangles) and Ising 
(squares) noise in an interacting system with a = 1, A = 
0.5 knT. The lines are guides to the eye. The size of this 
particular system was L = 32. 



Appendix C: Supporting simulations 
1. Polarized Ising randomness 

In order to confirm the picture that it is indeed the 
sites with Jp > that restrict the diffusion by acting as 
barriers to the anyons, we have determined the lifetime 
T of an interacting system (L = 50, a = Q, A = 0.5 ksT, 
a = 5 ksT) as a function of the Ising polarization P. The 
latter is defined as P = 1 — 2ri, where rj is the fraction of 
sites with negative onsite energy. 

Starting from P = —1, i.e., Jp = —a for all sites p, 
the lifetime moderately increases as more and more pos- 
itive sites are randomly added (increasing P). Around 
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FIG. 11. (color online) Lifetime as a function of Ising po- 
larization P. The parameters for these systems are L = 50, 
a = 0, A = 0.5 ksT, a = bknT. The dashed line is a guide 
to the eye. 



P = 0, where there is an equal number of sites with pos- 
itive and negative onsite energies, the hfetime drastically 
increases by about 2 orders of magnitude. At this point, 
large connected areas of sites with Jj, < can no longer 
exist, such that the anyons can move freely only within 
areas each consisting of just a few negative sites. Conse- 
quently, the diffusion is drastically reduced. If and how 
the polarization at this threshold is related to the site and 
bond percolation thresholds of the square lattice, which 
are rj 0.5927 and 0.5 (see, e.g., Ref. |25j), respectively, 
is not completely clear. 



2. Artificial cutoff of number of anyons 



anyons by simulating a non-interacting system with an 
artificial maximal number of anyons. This data is shown 
in Fig.[T2j Despite the absence of interaction, the lifetime 
of encoded states is still growing with increasing Ising 
disorder strength, hence clearly demonstrating that this 
effect is caused solely by the disorder. Furthermore, the 
lifetime saturates for large cr, because the energy barriers 
posed by the sites with Jp = -t-cr are essentially infinitely 
high and increasing them further bears no more advan- 
tage. The observed saturation also confirms that it is 
indeed the growing number of anyons that is responsible 
for the subsequent decrease in lifetime at large a in the 
data presented in the main text (where the number of 
anyons was not artificially restricted). 
12.5 




FIG. 12. (color online) Lifetime as a function of Ising disorder 
strength in a non-interacting system of size L = 32 with an 
artificially engineered maximal number of anyons equal to 20. 
The solid line is a guide to the eye. Inset: The error-corrected 
logical Z operator as a function of time for different o yielding 
the lifetimes shown in the main plot. 



We can support the claim that the only relevant effect 
of the repulsive interaction is to reduce the number of 
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